# source("./DATA_IN/source/llibreries_WGCNA.R")
source("/home/albert/albert/DATA_IN/source/llibreries_WGCNA.R")

## GRUP
grup<-"SA"
# Arrays
## 1. Lectura ja normalitzada

load("./DATA_IN/source/expression_set.RData")
dim_inicial<-dim(data_rma)

## 2. definciio fenotip
if(grup=="TUMOR"){
  g<-"T"
} else if(grup!="TUMOR"){
  g<-"S"
} else{g<-NA}
cond1<-data_rma@phenoData@data$tumor==g&data_rma@phenoData@data$position=="HIGH"
cond2<-data_rma@phenoData@data$tumor==g&data_rma@phenoData@data$position=="LOW"
data_rma$condicio<-rep("condicio",length(data_rma$index))
data_rma$condicio[cond1]<-paste0(grup,"_HIGH")
data_rma$condicio[cond2]<-paste0(grup,"_LOW")
data_rma_sa<-
  data_rma[,data_rma$condicio!="condicio"]

# Eliminar pacient 222
data_rma_sa<-data_rma_sa[,data_rma_sa@phenoData@data$sample!="322"]
datatable_mod(data_rma_sa@phenoData@data)

## 3. Filtrar per sd o per coevient de variació

library(genefilter)
dim(data_rma_sa)

ffun <- filterfun(cv(0.05,Inf))
t.fil <- genefilter(data_rma_sa,ffun)
# apply filter, and put expression back on log scale

dim( data_rma_sa[t.fil,])

data_rma_sa <- data_rma_sa[t.fil,]

# vars <- apply(data_rma_sa, 1,sd)
# sd_cut<-0.6
# print(paste0("Filtrat per sd= ",sd_cut))
# data_rma_sa<-(data_rma_sa[vars > quantile(vars, sd_cut), ] )
dim_1<-
dim(data_rma_sa)
print(paste0("Gens descartats= ",(dim_1-dim_inicial)[[1]]))

## Transformar a matriu d'expressió
data_rma_sa<-data_rma_sa[,order(data_rma_sa$position,decreasing = T)]
data_rma_sa_expr<-exprs(data_rma_sa)
datExpr0<-t(data_rma_sa_expr)
dim(datExpr0)

rownames(datExpr0)<-paste0(data_rma_sa@phenoData@data$sample,"_",data_rma_sa@phenoData@data$tumor,"_",data_rma_sa@phenoData@data$position)

# lipids
## Lectura

lipids<-read.csv("/home/albert/albert/APP_SA_srv_albert/WGCNA/DATA_IN/source/allTraits_new_v2_mod.csv")
lipids<-lipids %>% filter(delete!="si",tumor==g)
lipids$arrayID<-paste0(lipids$arrayID,"_",lipids$tumor,"_",lipids$pos)
datatable_mod(lipids)

groups = data_rma_sa@phenoData@data$condicio
groups<-as.factor(groups)
f = factor(groups)
design = model.matrix(~ 0 + f)
data.fit = lmFit(data_rma_sa_expr,design)

contrast.matrix = makeContrasts(paste(colnames(design)[2],"-",colnames(design)[1]),levels=design)
data.fit.con = contrasts.fit(data.fit,contrast.matrix)
data.fit.eb = eBayes(data.fit.con)
tab = topTable(data.fit.eb,coef=1,number=Inf,adjust.method="BH",genelist = data_rma_sa@featureData@data$SYMBOL)

write.csv(tab,"./RESULTATS/TAULES/FC.csv")

options(stringsAsFactors = FALSE)

gsg = goodSamplesGenes(datExpr0, verbose = 0);
if(gsg$allOK){print("Totes les mostres i gens passen el control de qualitat de WGCNA")}

if (!gsg$allOK){
# Optionally, print the gene and sample names that were removed:
if (sum(!gsg$goodGenes)>0)
printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")));
if (sum(!gsg$goodSamples)>0)
printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
# Remove the offending genes and samples from the data:
datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}
dim(datExpr0)

# Guardar la matriu d'expressió
save(datExpr0,file = "./RESULTATS/OBJECTES_R/datExpr0")

sampleTree = hclust(dist(datExpr0), method = "average");
# Plot the sample tree: Open a graphic output window of size 12 by 9 inches
# The user should change the dimensions if the window is too large or too small.
sizeGrWindow(12,9)

png("./RESULTATS/PLOTS/sampleClustering.png",width = 800,height = 800,res = 150)
par(cex = 0.6);par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
cex.axis = 1.5, cex.main = 2)
.close<-dev.off()
knitr::include_graphics("./RESULTATS/PLOTS/sampleClustering.png")

arraySamples = rownames(datExpr0)
arraySamples[!arraySamples%in%
lipids$arrayID]

traitRows = match(arraySamples, lipids$arrayID)

datTraits = lipids[traitRows, c(-3,-4,-5,-6,-7)]
rownames(datTraits)<-datTraits$arrayCel
datTraits$tumor<-gsub("S",0,datTraits$tumor)
datTraits$tumor<-gsub("T",1,datTraits$tumor)
datTraits$tumor<-as.numeric(datTraits$tumor)
datTraits$pos<-gsub("HIGH",1,datTraits$pos)
datTraits$pos<-gsub("LOW",0,datTraits$pos)
datTraits$pos<-as.numeric(datTraits$pos)


# Re-cluster samples
sampleTree2 = hclust(dist(datExpr0), method = "average")

# Convert traits to a color representation: white means low, red means high, grey means missing entry
traitColors = numbers2colors(datTraits, signed = F,lim = c(0,2));
# Plot the sample dendrogram and the colors underneath.
png("./RESULTATS/PLOTS/sampleClustering_2.png",width = 1000,height = 1000,res = 100)
plotDendroAndColors(sampleTree2, traitColors,
groupLabels = names(datTraits),
main = "Sample dendrogram and trait heatmap")
.close<-dev.off()
save(datExpr0, datTraits, file = paste0(grup,"_lipids-01-dataInput.RData"))
knitr::include_graphics("./RESULTATS/PLOTS/sampleClustering_2.png")

# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=50, by=2))
# Call the network topology analysis function


if(!file.exists("./RESULTATS/OBJECTES_R/pickSoftThreshold")){
sft = pickSoftThreshold(datExpr0,
                        powerVector = powers,
                        corFnc = "bicor",
                        networkType = "unsigned",
                        verbose = 0)
save(sft,file = "./RESULTATS/OBJECTES_R/pickSoftThreshold")
}else{load("./RESULTATS/OBJECTES_R/pickSoftThreshold")}

# Plot the results:

# Scale-free topology fit index as a function of the soft-thresholding power

png("./RESULTATS/PLOTS/beta.png",width = 1000,height = 1000,res = 150)
par(mfrow = c(1,2));
cex1 = 0.9;

plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
  # this line corresponds to using an R^2 cut-off of h
abline(h=0.9,col="red")
# # Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")

.close<-dev.off()

beta<-12
knitr::include_graphics("./RESULTATS/PLOTS/beta.png")

if(file.exists("./RESULTATS/OBJECTES_R/net")){
load("./RESULTATS/OBJECTES_R/net")  
}else{

net = blockwiseModules(datExpr0, power = beta, maxBlockSize = 20000,
                       TOMType = "unsigned",
                       corType = "bicor",
                       minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "./RESULTATS/OBJECTES_R/net", 
                       verbose = 3)

}

save(net,file="./RESULTATS/OBJECTES_R/net")

# open a graphics window
sizeGrWindow(12, 9)
# Convert labels to colors for plotting
mergedColors = labels2colors(net$colors)
# Plot the dendrogram and the module colors underneath
png("./RESULTATS/PLOTS/module_colors.png",height = 1000,width = 1000,res=150)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
.close<-dev.off()
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]]
save(MEs,net, moduleLabels, moduleColors, geneTree,file = "./RESULTATS/OBJECTES_R/net_auto")

knitr::include_graphics("./RESULTATS/PLOTS/module_colors.png")

# Define numbers of genes and samples
nGenes = ncol(datExpr0);
nSamples = nrow(datExpr0);
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datExpr0, moduleColors)$eigengenes
MEs = orderMEs(MEs0)

moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);


moduleTraitCor<-
moduleTraitCor[order(moduleTraitCor[,"pos"],decreasing = T),]

sizeGrWindow(10,6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(datTraits),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = greenWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = paste("Module-trait relationships"))


cor_value<-c()
p<-c()

for (i in 1:dim(moduleTraitCor)[2]){
  cor_value[i]<-stats::cor(moduleTraitCor[,"pos"],moduleTraitCor[,i])
}
# Ordenar el lipid segons la correalció amb la posició
moduleTraitCor_ord<-(moduleTraitCor[,order(cor_value)])


moduleTraitPvalue = corPvalueStudent(moduleTraitCor_ord, nSamples);


# sizeGrWindow(10,6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor_ord, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor_ord)


# sizeGrWindow(10,6)
# Will display correlations and their p-values

# par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor_ord,
xLabels = colnames(moduleTraitCor_ord),
yLabels = rownames(moduleTraitCor_ord),
ySymbols = rownames(moduleTraitCor_ord),
colorLabels = FALSE,
colors = greenWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = paste("Module-trait relationships"))

# library(WGCNA)
# Define variable weight containing the weight column of datTrait
for(i in 1:length(colnames(datTraits))){
weight = as.data.frame(datTraits[,i]);
if(var(weight)!=0){
names(weight) = colnames(datTraits)[i]
# names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr0, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
dim(geneModuleMembership)
table(rownames(geneModuleMembership)==
data_rma_sa@featureData@data$PROBEID)

names(geneModuleMembership) = paste("MM", modNames, sep="");
head(geneModuleMembership)
names(MMPvalue) = paste("p.MM", modNames, sep="");
geneModuleMembership$SYMBOL<-
  data_rma_sa@featureData@data$SYMBOL[match(rownames(geneModuleMembership),data_rma_sa@featureData@data$PROBEID)]

geneModuleMembership$PROBEID<-
  data_rma_sa@featureData@data$PROBEID[match(rownames(geneModuleMembership),data_rma_sa@featureData@data$PROBEID)]

geneTraitSignificance = as.data.frame(cor(datExpr0, weight, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(weight), sep="");
names(GSPvalue) = paste("p.GS.", names(weight), sep="")
f<-1
for(f in 1:length(modNames)){
  module<-modNames[f]
column = match(module, modNames);
moduleGenes = moduleColors==module;

par(mfrow = c(1,1));
png(paste0("./RESULTATS/PLOTS/MODULS/",colnames(datTraits)[i],"_",modNames[f],".png"))
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   xlab = paste("Module Membership in", module, "module"),
                   ylab = "Gene significance for body weight",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = "black")
                   
dev.off()                   
simbols<-geneModuleMembership[moduleGenes, "SYMBOL"]
gens_moduls<-                 
cbind(geneModuleMembership=geneModuleMembership[moduleGenes, column],
      geneTraitSignificance=abs(geneTraitSignificance[moduleGenes, 1]),
      SYMBOL=simbols)

write.csv(gens_moduls,paste0("./RESULTATS/TAULES/MODULS/",colnames(datTraits)[i],"_",modNames[f],".csv"))
}
}
}

knitr::include_graphics(paste0("./RESULTATS/PLOTS/MODULS/",colnames(datTraits)[i],"_",modNames[f],".png"))

dissTOM = 1-TOMsimilarityFromExpr(datExpr0, power = 6)

# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM = dissTOM^7;
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA;
# Call the plot function
sizeGrWindow(9,9)
png("./RESULTATS/PLOTS/heatmap.png",width = 1000,height = 1000)
TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
.as<-dev.off()
knitr::include_graphics("./RESULTATS/PLOTS/heatmap.png")

# Recalculate module eigengenes
MEs = moduleEigengenes(datExpr0, moduleColors)$eigengenes
# Isolate weight from the clinical traits
weight = as.data.frame(datTraits$PI_36.4);
names(weight) = "PI_36.4"

MET = orderMEs(cbind(MEs, weight))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5);
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
= 90)

sizeGrWindow(6,6);
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
plotHeatmaps = FALSE)
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
plotDendrograms = FALSE, xLabelsAngle = 90)